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Abstract 

This paper deals with numerical solutions to a partial differential equation of frac- 
tional order. Generally this type of equation describes a transition from anomalous 
diffusion to transport processes. From a phenomenological point of view, the equa- 
tion includes at least two fractional derivatives: spatial and temporal. In this paper 
we proposed a new numerical scheme for the spatial derivative, the so called Riesz- 
Feller operator. Moreover, using the finite difference method, we show how to employ 
this scheme in the numerical solution of fractional partial differential equations. In 
other words, we considered an initial-boundary value problem in one dimensional 
space. In the final part of this paper some numerical results and plots of simulations 
are shown as examples. 

Key words: Anomalous diffusion, Fractional calculus, Finite difference method, 
Riesz- Feller operator, Boundary conditions 



1 Introduction 



In the last years fractional derivatives have found numerous applications in many 
fields of physics, mathematics, mechanical engineering, biology, electrical engineer- 
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ing, control theory and finance [1,2,3,4,5]. One can find interesting properties and 
interpretations of fractional calculus in [4,5,6,7,8,9]. Fractional calculus in mathe- 
matics is a natural extension of integer-order calculus and gives a useful mathemat- 
ical tool for modelling many processes in nature. One of these processes, in which 
fractional derivatives have been successfully applied, is called diffusion. Phenomena 
which deviate from classical diffusion are described in many papers. The deviation 
is called anomalous diffusion. This type of diffusion is characterized by the nonlinear 
dependence of the mean square displacement x (t) of a diffusing particle over time t: 
(x 2 (t)) ~ /c 7 i 7 for < 7 < 2. This is the opposite of classical diffusion where 
the linear dependence ( x 2 (i)) ~ k±t occurs. Analysing changes in the parameter 7 
it may be said that transport phenomena in systems exhibiting sub-diffusion have 
< 7 < 1, and 1 < 7 < 2 in the systems exhibiting super-diffusion. However the 
dependence ( x 2 (t)} = 00 for t — > 00 characterizes rare but extremely large jumps of 
a diffusing particle - known as Levy motion or Levy flights [10]. Classical diffusion 
follows Gaussian statistics and Fick's second law for running processes at time t, 
whereas anomalous diffusion follows non-Gaussian statistics or can be interpreted 
as the Levy stable densities. The ever increasing amount of literature in which the 
behaviour of anomalous diffusion is observed presents some examples of this diffu- 
sion in: disordered vortex lattice for superconductors [11], supercooled liquids and 
glasses [12], disordered fractal media [13], liquid crystal polymers [14] and many 
others. 

Using equations with integer-order derivatives to model anomalous diffusion in the 
above processes may not reflect their real real behaviour Therefore anomalous diffu- 
sion is described by a spatio-temporal fractional partial differential equation in which 
classical spatial and temporal derivatives are replaced by derivatives of fractional or- 
der. The anomalous diffusion equation includes the classical diffusion equation, and 
therefore, this equation has a general form. Moreover, the anomalous diffusion equa- 
tion may also describe wave propagation or advection processes. 

Equations of anomalous diffusion with time and/or space fractional derivatives 
have been proposed and analysed by numerous authors, for example Nigmatullin 
in [15,16], Bouchaud in [17], Wyss and Schneider in [18,19,20], West in [21], Goren- 
flo, Fabritiis and Mainardi in [22], Gorenflo and Mainardi in [23], Mainardi in [24], 
Metzler and Klafter [25], Hilfer in [1] and recently by Agrawal in [26]. Nevertheless, 
the theoretical analysis and numerical methods applied to solve fractional diffusion 
equations present difficulties. In many papers, the autors have considered and solved 
problems in the infinite domain. Hilfer [1] and Klafter and Metzler [27] described the 
analytical solution to these equations in terms of Fox's ^/-function. In [26] Agrawal 
presented an analytical solution over time for the anomalous diffusion equation with 
boundary conditions of the first kind. He based his approach on the Laplace trans- 
form in terms of the Mittag-Leffer function. However, the numerical approximation 
for the series of expansions of these functions are a little problematic, especially for 
greater values of function arguments. 

Some numerical methods such as the finite difference method (FDM) and the fi- 
nite element method (FEM) are more suitable for solving the anomalous diffusion 
equation in more general (non-linear) form. The numerous works by Gorenflo and 



Mainardi [22,23], Podlubny [5] and many others should be noted. The difference 
scheme for fractional derivatives is based on the definition in the Grunwald-Letnikov 
form [5,9]. This from can make the scheme more flexible and straightforward. It is 
also a difficult task to solve the boundary value problem of these equations. Ciesielski 
and Leszczynski in [28] , Yuste in [29] proposed and analysed different cases of numer- 
ical solutions to the fractional diffusion equation where specific kinds of boundary 
conditions were selected. 

In this paper, we present a solution over space to a fractional partial differential 
equation in which the Riesz-Feller potential is included. In several papers [30,31] 
are shown some numerical schemes which didn't analyse originally proposed this 
potential. The authors introduced their own potentials which are suitable from a 
phenomenological point of view. Numerical treatments of space fractional differen- 
tial equations are not as popular in literature as the numerical solution to time frac- 
tional differential equations. This arises because the discretization of the Riesz-Feller 
potential is more difficult due to the appearance of singularities in this operator. On 
the basis of the Grunwald-Letnikov definition of fractional derivative, Gorenflo and 
Mainardi [22,23] proposed a method of discretization for this operator in the infinite 
domain. However, this descretization does not provide continuity for all values of 
the derivative order. Taking into account this inconvenience we shall propose an- 
other discrete form of the Riesz-Feller operator in order to restrict its continuity for 
arbitrary values of the derivative order. We use a new approach to solve numerically 
an equation of anomalous diffusion. 



2 Mathematical model 



In this paper we consider an equation in the following form 

-C(x,t) = K a ——C(x,t), t>0, xeR, (1) 
at o\x\ e 

where C(x, t) is a field variable, (d a /d \x\@ ) C(x, t) is the Riesz-Feller fractional oper- 
ator [32,25,9], a is the real order of this operator, K a is the coefficient of generalized 
(anomalous) diffusion with the unit of measurement m a /s. According to [22,23] the 
Riesz-Feller fractional operator for < a < 2, a ^ 1, and for the one-variable 
function u(x) is defined as 

= x DeU (x) = - [c L (a, 9) -ooDSu (x) + c R (a, 9) x D% 00 u {x)\ , (2) 

o \x\ e 
where 

[-ooJ2r a «(*)], (3) 

( i \ m 
-J [ x I™- a u(x)], (4) 



with m G N, m — 1 < a < m as the left-side and right-side Riemann-Liouville 
fractional derivatives. In formula (2) coefficients cl (ct,Q), cr {ol,9) (for < a < 2, 
and for \9\ < min (a, 2 — a)) have the following forms 

(a-e)7r 
sin 

c L (a,9) = — — - — , c R (a,0) 

sm (cnr) 

In expressions (3) and (4) the fractional operators -^l^u (x) and x I^u (x) are 
defined as the left- and right-side of Weyl fractional integrals [22,23,4,5,9]. Thus we 
have 

Assuming a = 2 in (1) we obtain classical diffusion, i.e. the heat transfer equation. 
On the other hand, the classical transport equation is obtained if a = 1 and the 
parameter of skewness 9 in (5) admits extreme values. Taking into consideration the 
above changes in parameter a we assume its variations within the range < a < 2. 
Analysing in (1) the behaviour of parameter a (a < 2) we observed some transition 
between the transport and propagation processes. 

We used the Green's function [23] in the analytical solution of Eq. (1). However, 
Eq. (1) needs a numerical solution when an additional non-linear term may occur. 
Some numerical methods used to solve fractional partial differential equations can be 
found in [22]. However, these methods apply the infinite domain without boundary 
conditions. 

In this work, we will consider Eq. (1) in the one dimensional domain £1 : L < x < R 
with boundary- value conditions of the first kind (Dirichlet conditions) as 

'x = L: C(L,t)=g L (t) 
< for t > 

x = R : C(R,t)=g R (t) 
and with the initial-value condition 

C(x,t)\ t=Q = c (x) 



3 Numerical method 

According to the FDM [33,34] we consider a discrete form of Eq. (1) both in time and 
space. In our previous work [28] we solved numerically an anomalous diffusion equa- 
tion similar to (1) where only the time- fractional derivative was taken into account. 



sm ■ 



{a + 9)iT 



sin (air) 



(5) 



(6) 
(7) 



(8) 



(9) 



We called this method the Fractional FDM (FFDM). Extending our considerations, 
we can say that the solution to equation (1) needs proper approximation of the 
Riesz-Feller derivative (2) in numerical schemes. 

Here, we introduce a definition of the fractional derivative in the Caputo form [35,5] 
as 

c a D«u{x) = a I™-«u^(x)i (10) 



c ;D%u{x)= x I™- a u^{x), (11) 

where m G N, m — 1 < a < m, a, 6 G R, a < x, b > x and are first and 

second derivatives, for m = 1,2. The following relations are present between the 
Riemann-Louville and the Caputo forms: 



m— 1 



d k 



a D-u{x) = c a D-u{x) + Y,^u{x 



k=0 
m—1 



d k 



^u{x)= c x D^u{x) + Y J ^ f«( 



k=0 



(x — a) k a 
x=a T(k-a+lY 

(b-x) k ~ a 
x=b T(k-a + iy 



(12) 
(13) 



Based on the assumptions |lim a ^ > _ 00 (d k /dx k ^ u (x) \ \ < oo and as well 
llim^oo (d k /dx k ) u(x)\ _,| < oo for k = 0, ...,m — 1, the terms occuring on the 
left sides of (12) and (13) tend to zero. Thus, when the lower/upper limit of inte- 
gration tends to minus/plus infinity we have 



- 00 D%u(x) = _£DZu(x) and x D° 00 u(x)=^D« 00 u(x). 



C r,a 



(14) 



3. 1 Approximation of the Riesz-Feller derivative 



As we start numerical analysis from the discretization of operators (10) and (11) 
respectively, we introduce a homogenous spatial grid — oo < . . . < Xj_2 < Xi-i < 
Xi < Xi + \ < Xi + 2 < . . . < oo with the step h = x^ — Xk-\- We denote the value of 
function u (x) at the point x^ as u& = u (x^), for fcsZ. We take into account only the 
function of one variable in order to simplify notations and we denote cl = cl (a, 9), 
C R = CR ( a i I n accordance with changes in parameter a in (1) we distinguish two 
cases of discrete approximation of the Riesz-Feller derivative. 

The first case includes changes in parameter a for the range < a < 1. We rewrite 
operator (2) using the Caputo definition as 



D$u{xi) 



C n « 

X 



CL 



r( 



I -a) I 



u'(0 



d£ - c R - 



(x,-O a s r(i 



oo 



r(l-a) 



-1 
r(i-a) 



OO ^i— * oo x i+k+l 



"'(0 



(15) 



x i — k 



Xi + k+l 



x i—k — l x i + k 



where u • and u^- are difference schemes which approximate the first derivative of 
integer order in the intervals [xj-±,Xj] and [xj,Xj+\], respectively. We propose the 
following weighed forms of these difference schemes as 



Ua = - 



Uj 



h 



Ai) (uj - Uj-i) + Ai (uj+i - Uj) 



— [Aiu j+ i + 2 (1 



Ai) Uj + (Ai 



h 

■ 2)v 



■j'-U 



(16) 



-/ 



1 / uj + i - uj (1 - Ai) (uj+i - Uj) + Ai (uj - Uj_i) 



= -L [(2 - Ax) u j+1 + 2 (Ai - 1) + (-Ax) u^x] , (17) 

where Ax = Ax (a, 6) = a — \6\, Ax € [0,1]. We introduce these formulae because 
we want to obtain various transitions between the difference schemes which are 
connected with the first derivative of integer order. For example, after putting Ax = 1 
into (16) and (17) we obtain wide known the central-difference approximation of 
first derivative, and after putting Ax = we get the backward- (16) or forward- (17) 
difference schemes. In this way we would like to avoid the problem of singularity 
in the operator (2) for a — > 1~ and 6 = (in further calculations of the discrete 
form of the Riesz-Feller operator we obtain finite values of the coefficents for this 
case) . Simultaneously, we can obtain classical schemes used for hyperbolic equations 
of first order for a — > 1~ and 6 — > 1^, e.g. (ui + \ — ui) /h and (uj — -u^x) /h. 

Denoting by 



k r(i-a) J (x t -0 a r(i-a) 7 (e-^) a ? 

— fc — 1 ^i + fc 

~ h U2^a) ' (18) 



we have 

XiDg Ui 



CL 12 2h ^ lUi - k+1 + 2 ( 1 ~ A u i~k + (Ai - 2) Uj_ fc _i] 4" 



(a) 



k=0 

oo 



CR 12^h t( 2 ~ Al ) + 2 (Ai - 1) ^i+fc + (-Ai) u i+k -i] 



.(«) 



fc=0 



Finally we can write the discrete form of (2) as 



fc=— oo 



where coefficients = u>fc (a, 0) have the following form: 



-1 



x < 



2r (2 - a) 

' ((|fc|+2) 1 - Q A 1 + (|A;| + l) 1 - Q (2-3A 1 ) 

+ {kl 1 - (3Ai - 4) + (|fc| - l) 1 "" (2 - Ai)) c L for A; < -2, 
(3 1 - Q Ai + 2 1 -" (2 - 3Ai) + 3Ai - 4) c L + Aic fi 
(2 1 ~ Q A 1 - 3Ai + 2) ( CL + c fl ) 
(3 1 - Q Ai + 2 1 -" (2 - 3Ai) + 3Ai - 4) c K + Aic L 
((fc + 2) 1 - a Ai + (A: + l) 1 - a (2-3Ai) 

+k 1 ~ a (3Ai - 4) + (jfe - l) 1 "" (2 - Ai)) c R for A; > 2. 



for A; = 


-1 


for k = 


o, 


for At = 


1, 



(19) 



(20) 



(21) 



The second case involves changes in parameter a for the range 1 < a < 2. As in 
previous calculations, we write operator (2) using the Caputo form as 



C x Wu{ Xi ) = - 



CL 



r(2 



Xi oo 

-a) J (Si-fl- 1 * +CR Y{2-a)J - x^'^ 



r(2 - a) 



iff i 



(22) 



where and Uj are difference schemes of the second derivative of integer order 
which we approximate by the following formulae 



If Uj + i-2v,j+Uj-i (l-\ 2 )(uj + i-2uj+Uj-i) + \ 2 {uj-2uj-i+Uj- 2 y' 



+ 



h 2 



2h 2 



[(2 - A 2 ) u j+1 + (3A 2 - 4) Uj + (2 - 3A 2 ) Uj_i + A 2 ^- 2 ] , 



(23) 



~" _ 1 ( Uj+i-2uj+Uj-i (l — \ 2 )(uj + i — 2uj+Uj-i) + \ 2 (uj +2 — 2uj + i+Uj)\ 
= ^2 [A2«j+2 + (2 - 3A 2 ) u J+1 + (3A 2 - 4) Uj + (2 - A 2 ) u^] , (24) 

where A 2 = A 2 (a, 9) = 2-(a + \6\), X 2 € [0, 1]. By putting A 2 = into (23) and (24) 
we determine the classical central-difference scheme, and for A 2 = 1 we obtain the 
backward/forward four-point discretizations of the second derivative of integer order. 



Denoting 



,(«) 



1 



x i — k 



1 



= h 



T(2-a) J (xi-0 

x i — k — l 



a—l s 



1 



1 



T(2-a) J (S-Xi) 

x i + k 



a—l s 



r (3 - a) 



(25) 



we have 



XiDgUi PS 



1 



fc=0 



X 1 

+A 2 ^_ fe _ 2 ] 4 a) + c R ^2 [ A 2^i+fc+2 + (2 - 3A 2 ) u i+k+1 



k=0 



+ (3A 2 - 4) u i+k + (2 - A 2 ) Uj+fe-i] ^ 



(a) 



(26) 



Finally we can write the discrete form of expression (2) as 



j oo 



(27) 



where coefficients Wk = w& (a, 0) are 



x \ 



-1 

2r (3 - a) 

' ((|fc| + 2) 2 ~ Q (2 - A 2 ) + (|fc| + lf- a (4A 2 - 6) + \k\ 2 ~ a (6 - 6A 2 ) 

+ (|fc| - l) 2 ~ a (4A 2 - 2) + (|fc| - 2) 2 " Q (-A 2 )) c L for k < -2, 

(3 2 " Q (2 - A 2 ) + 2 2 - a (4A 2 - 6) - 6A 2 + 6) c L + (2 - A 2 ) c R for fc = -1, 
( 2 2-a ( 2 _ A 2 ) + 4A 2 - 6) (c L + Cfl) for fc = 0, 



3 2-Q 



(2 - A 2 ) + 2 2 " a (4A 2 - 6) - 6A 2 + 6) c R + (2 - A 2 ) c L for fc = 1, 



((fc + 2) 2 - a (2 - A 2 ) + (k + 1) 2 ~ Q (4A 2 - 6) + k 2 - a (6 - 6A 2 ) 



+ (k - l) 2 - a (4A 2 - 2) + (k - 2) 2 ~ a (-A 2 )) c R 



for jfe > 2. 



(28) 



Summarising the above calculations, we presented difference schemes for the Riesz- 
Feller fractional derivative. It should be noted that expressions (20) and (27) are 
represented by the weighted sum over discrete values of function u at all the node's 
points. If index k tends to 0, i.e. in the nearest proximity of an arbitrary point Xj, 
one may observe higher values of Wk ■ Whereas values Wk decrease to zero for those 
nodes furthest from the point x«. Table 1 shows sample values of wj. calculated 
for different values of parameter a. Here we assumed 9 = 0, cl = cr, and some 
symmetry in coefficients Wk = w\k\, for k ^ 0. It should be noted that coefficients 
Wk(2,0) are identical as to those for the central difference approximation for the 
second derivative, for a = 2. Assuming the skewness parameter to be 9 = 1 ± and 
a = 1 ± we can see that coefficients Wfc(l ± ,l ± ) tend to values represented by the 
backward/forward difference scheme for the first derivative. 



Table 1 

Coefficients (a, 0) being dependent on parameter a and k respectively 



k 


a 




0.1 


0.5 


1- 


1.5 


2 





-0.993029 


-0.963132 


-0.857606 


-1.498970 


-2 


1 


0.041819 


0.170296 


0.253710 


0.574964 


1 


2 


0.022853 


0.067624 


0.064577 


0.125442 





3 


0.014264 


0.036213 


0.029047 


0.020048 





4 


0.010322 


0.023595 


0.016789 


0.009118 





5 


0.008054 


0.016974 


0.010996 


0.005125 





fO 


0.003751 


0.006116 


0.002926 


0.000906 






3.2 Complete fractional FDM 



Although discretization of the Riesz-Feller derivative in space was proposed, we 
describe the FDM for the equation of anomalous diffusion (1) in this subsection. 
Here, we restrict the numerical solution to only one-dimensional space. In comparison 
to the classical diffusion equation where discretization of the second derivative over 
space is approximated by the central difference scheme, the anomalous diffusion 
equation requires generalized schemes given by formulae (20) and (27), respectively. 
Boundary conditions have a direct influence to the numerical solution not only on 
boundary nodes but also in internal nodes of the domain. 

We introduce a temporal grid = t° < t l < . . . < v < t? +1 < . . . < oo with the 
grid step At = t^ +1 — P . At a point x k at the moment of time P we denote the 
function C (x, t) as C( = C (x k ,t f ) , for k G Z and / G N. 



3.2.1 Pure initial value problem 

In the explicit scheme of FDM we replaced (1) with the following formula 

c f+i _ c f 1 °° 

' = K ^ E c Lk w k- 



At 



(29) 



k= — oo 



After simplifications we obtained the final form as 



= E ^U> 



(30) 



k=— oo 



where coefficients pk = Pk (a, 0) are 



' At 

1 + K a - — wo for k = 0, 



Pk 



At 



h a 



(31) 



K a —w k for k ^ 0. 



Now we calculate the sum of all coefncent values p k : 



^ p k = l + K a — ^ Wk = l + K a — 

k=—oo k=— oo 



k=l 



(32) 



Substituting values from expressions (21) as well as (28) into w k and making many 
calculations in both cases we finally obtain 



E P k = 1 - 



(33) 



k=— oo 



In order to determine the stability of the explicit scheme (30) [36,37] the coefficient 
Po in (31) for k = should be positive (the other coefficients p k are non-negative for 



arbitrary values a and which is easy to prove) 

po = l + K a -^w > 0. (34) 

Thus we fixed the maximum length of the time step At, substituting into u>o ex- 
pressions (21) and (28) respectively, as 

T(2 -a) 



At < 



-h° 



2h c 



K a w K a (c L + c R ) 



2 1 -°Ai -3Ai + 2 
T (3 - a) 

I 2 2 -« (2- A 2 ) +4A 2 



6 



for < a < 1, 



for 1 < a < 2. 



(35) 



Moreover, the initial condition (9) is introduced directly to every grid node at the 
first time step t = t°. This determines the initial values of function C as 



(36) 



It is not easy to apply an implicit scheme in unbounded domains because it generates 
infinite dimensions for all matrices. Therefore, one usually seeks improved difference 
equations within an explicit scheme. 



3.2.2 The boundary-initial value problem 

The numerical solution (30), which included the unbounded domain — oo < x < oo, 
has no practical implementations in computer simulations. 

Here we try to solve this problem in the finite domain £1 : L < x < R with boundary 
conditions (8). We divide domain 0, into sub-domains with the step h = (R—L) /N. 
Figure 1 shows the modified spatial grid. Here, we can observe additional 'virtual' 

L R 



X-\ Xo X; Xn Xn+i-" 

Fig. 1. Grid nodes over space. 

points in the grid located outside the lower and upper limits of the domain Q. 
In order to introduce the Dirichlet boundary conditions, we propose a numerical 
treatment which assumes the same values of function C outside the domain limits 
as the values predicted on boundary nodes xq and xn- 



C (x k ,t) 



C(x ,t) = g L (t) for£;<0, 
C(x N ,t) = g R (t) iork>N. 



(37) 



Based on previous considerations we need to modify expressions (20) and (27) for 
the novel discretization of the Riesz-Feller derivative. Thus we have 



x .DgC (Xi, t) W — 



N- 



^2 C (x i+ k,t) w k + gi (t) s Li + g R (t) s RN . 



(38) 



for i = 1, . . . , N — 1, where 



SLj = SLj (a, #) = ^ ^fc 



= SRj (a, 0) = ^2 w k 
k=j+i 



cl ■ r j for < a < 1, 

C£ • rj for 1 < a < 2, 

cr - rj for < a < 1, 

cr -rj for 1 < a < 2 



(39) 
(40) 



and 



(j + 2f~ a Ai + (j + I) 1 -* (2 - 2Ai) + j l ~ a (Xi - 2) 



■l-a 



2T (2 - a) 



(j + 2f- a (2 - A 2 ) + (j + 1)"-" (3A 2 - 4) + j (2 - 3A 2 ) + (j - 1)^° A 

2T (3 - q) 



(41) 



(42) 



Modyfing expression (29) by substituting (38) we obtain a finite difference scheme 
which is dependent on the weight factor cr. Here we assume 



9l +1/2 = 9 L (tf+^)= 9L (At (/+!)), 
^ +1/2 =^(* /+1/2 )=^( A *(/ + ^)) 



(43) 
(44) 



in order to simplify the numerical scheme. For internal nodes Xj, i = 1, . . . , N — 1 
we have 



At 



JV-i 



,k=—i 



/+1/2 /+1/2 

+ % S U + 9 R s RN _ 



and for boundary nodes xq, xn we denote 



U N - 9r 



(45) 



(46) 
(47) 



The method is explicit for cr = 1, partially implicit for < a < 1 and fully implicit 
for cr = 0. 



The above scheme described by expressions (45)- (47) can written in the matrix form 



as 



A C f+1 = B, 



(48) 



where 





1 








. 


. 














a_i 


1 + a 


ai 


a 2 . 


• (IN-3 


ajv-2 


ajv-i 




6i 




a_ 2 


a_i 


1 + a 


ai 


. ajv-4 


OAT-3 


OAf-2 




6 2 




a_ 3 


a_ 2 


a_i 


1 + a . 


• Oat_5 


aAr„4 


OAf-2 




^3 


A = 


a_4 


a_ 3 


a-2 


a_i . 




ajv-3 


OjV— 4 


, B = 


5 4 




a-Af+2 


O-Af+3 


a-JV+4 


O-JV+5 • 


• 1 + a 


ai 


a 2 




bN-2 




a-AT+i 




a-N+3 


a_Ar + 4 . 


. a_i 


1 + ao 


ai 




&JV-1 













. 


. 





1 




s £ +i/2 




















(4 



with 



At 

Oj = (a - 1) K Q - Wj for j = -iV + 1, . . . , iV - 1, 



5l + 9 R 2 srn-j + o- 2^ Q+fc^ 

fc=-j 

for j = 1,..., iV- 1. 



and C f+1 is the vector of unknown function values at time t^ +1 . 



(50) 



(51) 



A particular case of the scheme (45) is the explicit scheme for a = 1. This case may 
be simplified to 



a 



Tl — I 



for i = 0, 



C/ +fe Pfe + A" a ^ (5{ +1/2 SLi + 9 f R 1/2 s RN -^ for i = 1, . 

for i = AT, 



(52) 

where is defined by formula (31). We can observe that boundary conditions in- 
fluence on all the values of the function at every node. Unlike the classical second 
derivative which is approximated locally, the Riesz- Feller and other fractional deriva- 
tives accumulate all values of the function at the domain points. 

The skewness parameter has a significant influence on the solution. Assuming 
6 -► ±1± and a -► 1± one can obtain the classical hyperbolic equation called the first 
order wave equation (the transport equation). In this case our scheme tends to the 
known Euler's forward-time and backward-space/forward-space scheme. For a = 2, 
9 = our scheme is the same as the forward-time and central-space scheme [33,34]. 



4 Simulation results and their analysis 



In this section we present the results of calculation obtained by our numerical ap- 
proach. We try to simulate the evolutions of the probability density function over 
time for a £ (0.1,0.5, 1~, 1.5,2). Assuming the domain [—10,10] we divided this 
domain into 1000 subintervals (h = 0.02). We assumed the boundary conditions as 
gL (t) = gR (t) = 0. The initial condition u(x, + ) = 5(x) is approximated by C® 00 = 
1/h and Cf = 0, for i / 500. Fig. 2 shows graphically the probability density func- 
tion G a (x, t) over space after t = Is in the limited interval [—7, 7] (in the logarithmic 
scale). It should be noted that for a = 2 our solution roughly estimates the Gaussian 
probability density function: G2(x,l) = (2y/ir) 1 exp (— x 2 ). For a = 1^ this solu- 
tion becomes the Cauchy probability density function: G\(x,l) = 1/ (l + x 2 )). 




Fig. 2. Probability density functions over space for different values of parameter a. 



In the next two examples we show the influence of parameters a and 9 on the 
solution. In both examples we consider the domain : [0, 1] with initial condition 



c (x) 



10 for x £ [0.4, 0.6] , 

for x £ [0, 0.4) U (0.6, 1] 



(53) 



Fig. 3 shows the solution to function C (x, t) for a = 0.9 and 9 = —0.7 with boundary 
conditions gL (t) = ga (i) = on Fig. 3(a) and gL (t) = 10, gn (t) = on the 
Fig. 3(b), at different moments of time. Whereas Fig. 4 presents plots of C (x,t), for 
a = 1.6, 9 = -0.4 and g L (t) = 10, g R (t) = 0. 




In summary we proposed the FFDM for fractional diffusion equations in which the 
Riesz-Feller fractional derivative is included. We analysed the anomalous diffusion 
equation in linear form in order to compare numerical results with the analytical 
solution. We obtained implicit and explicit FDM schemes which may generalise 
classical schemes of FDM. Moreover, our solution for a = 2 is the same as the 
classical finite difference method. We hope that this numerical approach will be 
successfully applied to fractional partial differential equations having more complex 
forms, i.e. non-linear forms. 

Analysing the graphs included in this work, we observe that for the case a < 2 
(the Levy flight) diffusion is slower than classical diffusion (Brownian motion) at 
the initial time steps. Nevertheless, the probability density function generates a long 



0,0 0,2 0,4 x °' 6 °' 8 1 -° 

Fig. 4. Numerical solution to Eq. (1) over space for a = 1.6, 9 = —0.4. 

tail of distribution in the long time limit. This can be associated with rare and 
extreme events which are characterized by very large arbitrary values of particle 
jumps. 

Analysing changes in the skewness parameter 6 we observed interesting behaviour 
in the numerical solution to Eq. (1). For a — > 1 + and 9 — > ±1 + we obtained the 
first order wave equation. Assuming € (0, 1) (with restrictions to the order a) we 
generated a class of non-symmetric probability density functions. 

The proposed numerical scheme creates a bridge between Gaussian and Cauchy 
processes. Our scheme is also a bridge between diffusion and transport phenomena. 
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